/**
   Copyright 2008 Kirill Ignatiev

    This file is part of kerr-wavefronts.

    kerr-wavefronts is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    kerr-wavefronts is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with kerr-wavefronts.  If not, see <http://www.gnu.org/licenses/>.
 */
#include <cstdio>
#include <cmath>
#include <gsl/gsl_errno.h>
#include "solvegeods.h"
using namespace std;

#include "geodesic.h"

#define loop(i,n) for (int i = 0; i < n; ++i)

void geodesic::evolve_ode(double tmax) {
  // Evolve
  double h = 1e-3;
  int status = GSL_SUCCESS;
  while (status == GSL_SUCCESS && t < tmax) {
    status = gsl_odeiv_evolve_apply(ev, ctl, step, &sys, &t, tmax, &h, data);
  }
  if (status != GSL_SUCCESS)
    assert(false);
}

schw_geodesic::schw_geodesic(double r0, double phi0, double rdot0,
                             double phidot0) 
  : geodesic(r0,phi0) {
  l = phidot0 * r0 * r0;
  sys.params = &l;
  sys.function = &schw_f;
  //fprintf(stderr, "Schw black hole: rdot %f, phidot %f, beta %f, horizon %f.\n", rdot0, phidot0, l, 1.0); fflush(stderr);
  rdot0 = rdot0 / sqrt(1-1/r0);
  data[0] = schw_x_of_r(r0);
  data[1] = phi0;
  data[2] = rdot0;
  r = r0;
  t = 0;
  phi = phi0;
}

void schw_geodesic::evolve(double tmax) {
  if (r < 1+1e-3) return;
  evolve_ode(tmax);
  r = schw_r_of_x(data[0]);
  phi = data[1];
}

kerr_geodesic::kerr_geodesic(const kerr_param& p, double r0,
                             double phi0, double pr0)
  : geodesic(r0,phi0), param(p) {
  sys.function = &kerr_f;
  sys.params = &param;
  data[0] = kerr_x_of_r(param, r0);
  data[1] = phi0;
  data[2] = pr0;
  r = r0;
  t = 0;
  phi = phi0;
}
void kerr_geodesic::evolve(double tmax) {
  if (!not_on_horizon()) return;
  // Evolve
  evolve_ode(tmax);
  r = kerr_r_of_x(param, data[0]);
  phi = data[1];
}
